Risk Groups in Penile Cancer
Aim of the study
The aim of this study is to evaluate the accuracy of previously published risk groups systems for predicting inguinal nodal metastases in patients with penile carcinoma. Cases of invasive penile squamous cell carcinomas were stratified using the following systems: Solsona et al (J Urol 2001;165:1506), Hungerhuber et al (Urology 2006;68:621), and the proposed by the European Association of Urology (Eur Urol 2004;46:1), with low, intermediate, and high risk categories in each one of them. Metastatic rates and cancer-specific survival rates in our dataset were compared with previously reported results. Receiver-operator characteristic (ROC) analysis was carried out to compare accuracy in predicting final nodal status. We found that these risk groups systems may be useful for patients with low-grade superficial tumors and less accurate for evaluating patients with high-grade locally-advanced penile carcinomas. The results of this study may be useful for therapeutic planning of patients with penile squamous cell carcinomas.
Description of the repository
This repository contains the full statistical analysis of the dataset that was used for the article “Risk Groups Systems for Penile Cancer Management: A Study of 203 Patients with Invasive Squamous Cell Carcinoma”. The article has been accepted for publication in Urology. This repository also contains the following files:
- The final (approved) PDF version of the article, as submitted for publication
- The BibTeX file containing all the references cited in the article
- The R script that was used for analyzing the dataset and write the article, as well as the R scripts containing the functions written for plotting figures and tables
- The R Markdown file used for this report
- The figures included in this repository in PNG format
Data were analyzed using R version 3.1.1 “Sock it to Me” (R Foundation for Statistical Computing, Vienna, Austria). Results were written using RMarkDown in RStudio version 0.98.1102 and the knitr package version 1.9 by Yihui Xie.
Building the dataset for analysis
First we loaded the full dataset including 333 patiens with invasive penile squamous cell carcinoma.
Data <- read.csv("Article/PenisSCC_333.csv")
The full dataset is available at http://dx.doi.org/10.6084/m9.figshare.1290997, which also contains the dataset’s codebook. From the 333 patients we selected only those with total/partial penectomy in Procedure, excluded cases with missing values in the variables of Grade, pT and cN, and excluded patients who were lost at follow up in Outcome.
Data <- subset(Data, Procedure == "Total penectomy" | Procedure == "Partial penectomy")
Data <- Data[complete.cases(Data$Grade), ]
Data <- Data[complete.cases(Data$pT), ]
Data <- Data[complete.cases(Data$cN), ]
Data <- subset(Data, Outcome != "Lost at follow-up")
Data <- droplevels(Data)
We then recoded 3 variables, creating new variables for cN positivity (cN_Positive), cancer-related death (DOD), and final nodal status (Final_Nodal). Positive cases for Final_Nodalwhere those with lymph node metastasis in the groin dissection, local relapse during follow-up, or unfavorable outcome, including alive with disease and death by cancer.
Data$cN_Positive <- ifelse(Data$cN == "cN0", c("No"), c("Yes"))
Data$cN_Positive <- as.factor(Data$cN_Positive)
Data$DOD <- ifelse(Data$Outcome == "Died of cancer", c("Yes"), c("No"))
Data$DOD <- as.factor(Data$DOD)
Data$Final_Nodal <- ifelse(Data$Mets == "Yes", c("Positive"), c("Negative"))
Data$Final_Nodal[Data$Local == "Yes"] <- "Positive"
Data$Final_Nodal[Data$Outcome == "Alive with disease"] <- "Positive"
Data$Final_Nodal[Data$Outcome == "Died of cancer"] <- "Positive"
Data$Final_Nodal <- as.factor(Data$Final_Nodal)
We finally created the risk groups for the Solsona et al, EUA, and Hungerhuber et al systems, using the criteria provided in the reports.
# SOLSONA risk groups
Data$Solsona[Data$Grade == "Grade 1" & Data$pT == "T1"] <- "Low risk"
Data$Solsona[Data$Grade == "Grade 2" & Data$pT == "T1"] <- "Intermediate risk"
Data$Solsona[Data$Grade == "Grade 3" & Data$pT == "T1"] <- "Intermediate risk"
Data$Solsona[Data$Grade == "Grade 1" & Data$pT == "T2"] <- "Intermediate risk"
Data$Solsona[Data$Grade == "Grade 1" & Data$pT == "T3"] <- "Intermediate risk"
Data$Solsona[Data$Grade == "Grade 2" & Data$pT == "T2"] <- "High risk"
Data$Solsona[Data$Grade == "Grade 3" & Data$pT == "T2"] <- "High risk"
Data$Solsona[Data$Grade == "Grade 2" & Data$pT == "T3"] <- "High risk"
Data$Solsona[Data$Grade == "Grade 3" & Data$pT == "T3"] <- "High risk"
Data$Solsona <- factor(Data$Solsona, levels = c("Low risk", "Intermediate risk", "High risk"), ordered = TRUE)
# EAU risk groups
Data$EAU[Data$Grade == "Grade 1" & Data$pT == "T1"] <- "Low risk"
Data$EAU[Data$Grade == "Grade 2" & Data$pT == "T1"] <- "Intermediate risk"
Data$EAU[Data$Grade == "Grade 3" & Data$pT == "T1"] <- "High risk"
Data$EAU[Data$pT == "T2"] <- "High risk"
Data$EAU[Data$pT == "T3"] <- "High risk"
Data$EAU <- factor(Data$EAU, levels = c("Low risk", "Intermediate risk", "High risk"), ordered = TRUE)
# HUNGERHUBER risk groups
Data$Hungerhuber[Data$Grade == "Grade 1" & Data$pT == "T1"] <- "Low risk"
Data$Hungerhuber[Data$Grade == "Grade 2" & Data$pT == "T1"] <- "Low risk"
Data$Hungerhuber[Data$Grade == "Grade 1" & Data$pT == "T2"] <- "Intermediate risk"
Data$Hungerhuber[Data$Grade == "Grade 1" & Data$pT == "T3"] <- "Intermediate risk"
Data$Hungerhuber[Data$Grade == "Grade 2" & Data$pT == "T2"] <- "Intermediate risk"
Data$Hungerhuber[Data$Grade == "Grade 2" & Data$pT == "T3"] <- "Intermediate risk"
Data$Hungerhuber[Data$Grade == "Grade 3"] <- "High risk"
Data$Hungerhuber <- factor(Data$Hungerhuber, levels = c("Low risk", "Intermediate risk", "High risk"), ordered = TRUE)
Data analysis was carried out on this dataset, using the simpleR package.
library(simpleR)
Methodology
Data analysis is divided in 4 sections, as it follows:
Descriptive Statistics. All the variables included in the dataset were analyzed using bar plots, histograms, box plots, and one-way tables. Factor variables were described using absolute and relative percentages. Numeric variables were described using mean, standard deviation, median, interquartile range, minimum and maximum value.
Inferential Statistics: Statistical tests (Fisher’s exact test for categorical variables, Kruskal-Wallis test for numerical variables) were used to evaluate the association between Final Nodal Status and Cancer-Related Death as predictors and all the variables included in the dataset. A 2-tailed P value was reported in all instances. Reported statistics included absolute and relative percentages for categorical variables; and mean, standard deviation, median, interquartile range, minimum and maximum value for numeric variables, by final nodal status and by cancer-related death.
Survival Analysis. For all variables in the dataset survival curves by Final Nodal Status and by Cancer-Related Death were built using the Kaplan-Meier method and compared using the Mantel-Cox (log-rank) test. Numerical variables were splitted in 2 levels using the median as the cutoff point. A 2-tailed P value was reported in all instances.
Receiver-Operator Characteristic (ROC) Curve Analysis. ROC curves for predicting Final Nodal Status and Cancer-Related Death were plotted for all risk groups systems. ROC curves were compared against the baseline using the Mann-Whitney U test. Reported statistics included the area under the ROC curve (AUC) with tis 95% confidence interval and the 2-tailed P value from the Mann-Whitney test.
Here it follows the description of all the variables included in the analyzed dataset.
Surgical procedure for primary treatment
Var <- Data$Procedure
categorical.plot(Var)

descriptive.categorical(Var)
| Partial penectomy |
124 |
61 |
| Total penectomy |
79 |
39 |
Number of missing cases: 0 cases.
Histologic grade
Var <- Data$Grade
categorical.plot(Var)

descriptive.categorical(Var)
| Grade 1 |
52 |
26 |
| Grade 2 |
69 |
34 |
| Grade 3 |
82 |
40 |
Number of missing cases: 0 cases.
Histologic subtype
Var <- Data$Subtype
categorical.plot(Var, align = "h", left =8)

descriptive.categorical(Var)
| Adenosquamous |
3 |
1.48 |
| Basaloid |
9 |
4.43 |
| Mixed |
21 |
10.34 |
| Papillary |
6 |
2.96 |
| Sarcomatoid |
2 |
0.99 |
| Usual |
132 |
65.02 |
| Verrucous |
17 |
8.37 |
| Warty |
13 |
6.40 |
Number of missing cases: 0 cases.
Tumor located in glans
Var <- Data$Glans
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 3 cases.
Tumor located in coronal sulcus
Var <- Data$Sulcus
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 3 cases.
Tumor located in Foreskin
Var <- Data$Foreskin
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 3 cases.
Anatomical level of maximum tumor invasion
Var <- Data$Level
categorical.plot(Var, align = "h", left = 9)

descriptive.categorical(Var)
| Corpus cavernosum |
122 |
60.1 |
| Corpus spongiosum |
65 |
32.0 |
| Lamina propria |
9 |
4.4 |
| Preputial dartos |
7 |
3.4 |
Number of missing cases: 0 cases.
Tumor thickness in mm
Var <- Data$Thickness
numerical.plot(Var, label = "Tumor Thickness, mm")


descriptive.numerical(Var)
| Mean |
7.5 |
| Standard Deviation |
3.7 |
| Median |
7 |
| Interquartile Range |
5 |
| Mininum |
1 |
| Maximum |
20 |
Number of missing cases: 11 cases.
Tumor size in cm
Var <- Data$Size
numerical.plot(Var, label = "Tumor Size, cm")


descriptive.numerical(Var)
| Mean |
4.7 |
| Standard Deviation |
2.5 |
| Median |
4 |
| Interquartile Range |
2.8 |
| Mininum |
1 |
| Maximum |
16 |
Number of missing cases: 91 cases.
Patient’s age in years
Var <- Data$Age
numerical.plot(Var, label = "Patient's Age, years")


descriptive.numerical(Var)
| Mean |
53.9 |
| Standard Deviation |
14.1 |
| Median |
53 |
| Interquartile Range |
21 |
| Mininum |
24 |
| Maximum |
88 |
Number of missing cases: 0 cases.
Bilateral inguinal lymphadenectomy
Var <- Data$Lymphadenectomy
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 0 cases.
Tumor invasion of penile urethra
Var <- Data$Urethra
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 0 cases.
Vascular invasion in primary tumor
Var <- Data$Vascular
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 7 cases.
Perineural invasion in primary tumor
Var <- Data$Perineural
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 6 cases.
Pathological T stage
Var <- Data$pT
categorical.plot(Var)

descriptive.categorical(Var)
| T1 |
8 |
3.9 |
| T2 |
107 |
52.7 |
| T3 |
88 |
43.3 |
Number of missing cases: 0 cases.
Any tumor relapse (local, regional or systemic)
Var <- Data$Relapse
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 8 cases.
Local tumor relapse
Var <- Data$Local
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 158 cases.
Regional tumor relapse
Var <- Data$Regional
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 158 cases.
Systemic tumor relapse
Var <- Data$Systemic
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 158 cases.
Follow-up length in months
Var <- Data$FollowUp
numerical.plot(Var, label = "Follow-Up, Months")


descriptive.numerical(Var)
| Mean |
108.1 |
| Standard Deviation |
94.7 |
| Median |
86.8 |
| Interquartile Range |
152.4 |
| Mininum |
0.8 |
| Maximum |
433.9 |
Number of missing cases: 0 cases.
Patient’s outcome
Var <- Data$Outcome
categorical.plot(Var)

descriptive.categorical(Var)
| Alive with disease |
14 |
6.9 |
| Alive without disease |
89 |
43.8 |
| Died of cancer |
34 |
16.7 |
| Died of other causes |
66 |
32.5 |
Number of missing cases: 0 cases.
Clinical N stage
Var <- Data$cN
categorical.plot(Var)

descriptive.categorical(Var)
| cN0 |
108 |
53.2 |
| cN1 |
31 |
15.3 |
| cN2 |
59 |
29.1 |
| cN3 |
5 |
2.5 |
Number of missing cases: 0 cases.
Clinically positive inguinal lymph nodes
Var <- Data$cN_Positive
categorical.plot(Var)

descriptive.categorical(Var)
Number of missing cases: 0 cases.
Final nodal status
Var <- Data$Final_Nodal
categorical.plot(Var)

descriptive.categorical(Var)
| Negative |
132 |
65 |
| Positive |
71 |
35 |
Number of missing cases: 0 cases.
Solsona et al risk groups
Var <- Data$Solsona
categorical.plot(Var)

descriptive.categorical(Var)
| Low risk |
5 |
2.5 |
| Intermediate risk |
50 |
24.6 |
| High risk |
148 |
72.9 |
Number of missing cases: 0 cases.
European Association of Urology risk groups
Var <- Data$EAU
categorical.plot(Var)

descriptive.categorical(Var)
| Low risk |
5 |
2.46 |
| Intermediate risk |
2 |
0.99 |
| High risk |
196 |
96.55 |
Number of missing cases: 0 cases.
Hungerhuber et al risk groups
Var <- Data$Hungerhuber
categorical.plot(Var)

descriptive.categorical(Var)
| Low risk |
7 |
3.4 |
| Intermediate risk |
114 |
56.2 |
| High risk |
82 |
40.4 |
Number of missing cases: 0 cases.
Var2 <- Data$Final_Nodal
Final nodal status and surgical procedure
Var1 <- Data$Procedure
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| Partial penectomy |
95 |
72 |
29 |
40.8 |
| Total penectomy |
37 |
28 |
42 |
59.2 |
Final nodal status and histologic grade
Var1 <- Data$Grade
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| Grade 1 |
48 |
36.4 |
4 |
5.6 |
| Grade 2 |
46 |
34.8 |
23 |
32.4 |
| Grade 3 |
38 |
28.8 |
44 |
62.0 |
Final nodal status and glans location
Var1 <- Data$Glans
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
2 |
1.6 |
0 |
0 |
| Yes |
127 |
98.4 |
71 |
100 |
Final nodal status and coronal sulcus location
Var1 <- Data$Sulcus
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
95 |
73.6 |
51 |
71.8 |
| Yes |
34 |
26.4 |
20 |
28.2 |
Final nodal status and foreskin location
Var1 <- Data$Glans
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
2 |
1.6 |
0 |
0 |
| Yes |
127 |
98.4 |
71 |
100 |
Final nodal status and anatomical level
Var1 <- Data$Level
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| Corpus cavernosum |
68 |
51.5 |
54 |
76.1 |
| Corpus spongiosum |
50 |
37.9 |
15 |
21.1 |
| Lamina propria |
9 |
6.8 |
0 |
0.0 |
| Preputial dartos |
5 |
3.8 |
2 |
2.8 |
Final nodal status and tumor thickness
Var1 <- Data$Thickness
numerical.group.plot(Var1, Var2)

descriptive.numerical.group(Var1, Var2)
| Mean |
6.8 |
8.7 |
| Standard Deviation |
3.3 |
4.1 |
| Median |
7.0 |
8.0 |
| Interquartile Range |
3.0 |
6.0 |
| Minimum |
1.0 |
2.0 |
| Maximum |
20.0 |
20.0 |
Final nodal status and tumor size
Var1 <- Data$Size
numerical.group.plot(Var1, Var2)

descriptive.numerical.group(Var1, Var2)
| Mean |
4.4 |
5.2 |
| Standard Deviation |
2.4 |
2.7 |
| Median |
4.0 |
5.0 |
| Interquartile Range |
2.5 |
3.5 |
| Minimum |
1.5 |
1.0 |
| Maximum |
16.0 |
12.0 |
Final nodal status and patient’s age
Var1 <- Data$Age
numerical.group.plot(Var1, Var2)

descriptive.numerical.group(Var1, Var2)
| Mean |
54.3 |
53.0 |
| Standard Deviation |
14.5 |
13.4 |
| Median |
53.0 |
53.0 |
| Interquartile Range |
22.0 |
19.0 |
| Minimum |
24.0 |
27.0 |
| Maximum |
88.0 |
82.0 |
Final nodal status and urethral invasion
Var1 <- Data$Urethra
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
84 |
63.6 |
31 |
43.7 |
| Yes |
48 |
36.4 |
40 |
56.3 |
Final nodal status and vascular invasion
Var1 <- Data$Vascular
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
104 |
81.2 |
37 |
54.4 |
| Yes |
24 |
18.8 |
31 |
45.6 |
Final nodal status and perineural invasion
Var1 <- Data$Perineural
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
109 |
85.2 |
25 |
36.2 |
| Yes |
19 |
14.8 |
44 |
63.8 |
Final nodal status and clinical N stage
Var1 <- Data$cN
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| cN0 |
92 |
69.7 |
16 |
22.5 |
| cN1 |
13 |
9.8 |
18 |
25.4 |
| cN2 |
25 |
18.9 |
34 |
47.9 |
| cN3 |
2 |
1.5 |
3 |
4.2 |
Final nodal status and pathological T stage
Var1 <- Data$pT
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| T1 |
8 |
6.1 |
0 |
0.0 |
| T2 |
76 |
57.6 |
31 |
43.7 |
| T3 |
48 |
36.4 |
40 |
56.3 |
Final nodal status and (any) tumor relapse
Var1 <- Data$Relapse
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
126 |
95.5 |
24 |
38.1 |
| Yes |
6 |
4.5 |
39 |
61.9 |
Final nodal status and local relapse
Var1 <- Data$Local
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
6 |
100 |
31 |
79.5 |
| Yes |
0 |
0 |
8 |
20.5 |
Final nodal status and regional relapse
Var1 <- Data$Regional
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
0 |
0 |
9 |
23.1 |
| Yes |
6 |
100 |
30 |
76.9 |
Final nodal status and systemic relapse
Var1 <- Data$Systemic
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| No |
6 |
100 |
27 |
69.2 |
| Yes |
0 |
0 |
12 |
30.8 |
Final nodal status and Solsona et al risk groups
Var1 <- Data$Solsona
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| Low risk |
5 |
3.8 |
0 |
0.0 |
| Intermediate risk |
46 |
34.8 |
4 |
5.6 |
| High risk |
81 |
61.4 |
67 |
94.4 |
Final nodal status and EAU risk groups
Var1 <- Data$EAU
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| Low risk |
5 |
3.8 |
0 |
0 |
| Intermediate risk |
2 |
1.5 |
0 |
0 |
| High risk |
125 |
94.7 |
71 |
100 |
Final nodal status and Hungerhuber et al risk groups
Var1 <- Data$Hungerhuber
categorical.group.plot(Var1, Var2)

descriptive.categorical.group(Var1, Var2)
| Low risk |
7 |
5.3 |
0 |
0 |
| Intermediate risk |
87 |
65.9 |
27 |
38 |
| High risk |
38 |
28.8 |
44 |
62 |
# Defining outcome variable
Status <- Data$Final_Nodal
# Creating dicotomic variables from numerical variables for plotting
Thickness_Median <- factor(ifelse(Data$Thickness > median(Data$Thickness, na.rm = TRUE), c("Above Median Thickness"), c("Below Median Thickness")))
Age_Median <- factor(ifelse(Data$Age > median(Data$Age, na.rm = TRUE), c("Above Median Age"), c("Below Median Age")))
# By surgical procedure
with(Data, survival.plot(Procedure, FollowUp, Status, title = "Final Nodal Status by Surgical Procedure"))

# By histologic grade
with(Data, survival.plot(Grade, FollowUp, Status, title = "Final Nodal Status by Histologic Grade", position = "bottomright"))

# By glans invasion
with(Data, survival.plot(Glans, FollowUp, Status, title = "Final Nodal Status by Glans Invasion"))

# By coronal sulcus invasion
with(Data, survival.plot(Sulcus, FollowUp, Status, title = "Final Nodal Status by Coronal Sulcus Invasion"))

# By foreskin invasion
with(Data, survival.plot(Foreskin, FollowUp, Status, title = "Final Nodal Status by Foreskin Invasion"))

# By anatomical level
with(Data, survival.plot(Level, FollowUp, Status, title = "Final Nodal Status by Anatomical Level"))

# By median tumor thickness
with(Data, survival.plot(Thickness_Median, FollowUp, Status, title = "Final Nodal Status by Median Tumor Thickness"))

# By median patient's age
with(Data, survival.plot(Age_Median, FollowUp, Status, title = "Final Nodal Status by Median Patient's Age"))

# By urethral invasion
with(Data, survival.plot(Urethra, FollowUp, Status, title = "Final Nodal Status by Urethral Invasion"))

# By vascular invasion
with(Data, survival.plot(Vascular, FollowUp, Status, title = "Final Nodal Status by Vascular Invasion"))

# By perineural invasion
with(Data, survival.plot(Perineural, FollowUp, Status, title = "Final Nodal Status by Perineural Invasion"))

# By pathological T stage
with(Data, survival.plot(pT, FollowUp, Status, title = "Final Nodal Status by Pathological T Stage"))

# By clinical N stage positivity
with(Data, survival.plot(cN_Positive, FollowUp, Status, title = "Final Nodal Status by Clinical N Stage Positivity"))

# By Solsona et al risk groups
with(Data, survival.plot(Solsona, FollowUp, Status, title = "Final Nodal Status by Solsona et al Risk Groups", position = "bottomright"))

# By EAU risk groups
with(Data, survival.plot(EAU, FollowUp, Status, title = "Final Nodal Status by EAU Risk Groups"))

# By Hungerhuber et al risk groups
with(Data, survival.plot(Hungerhuber, FollowUp, Status, title = "Final Nodal Status by Hungerhuber et al Risk Groups"))

# Defining outcome variable
Status <- Data$DOD
# By surgical procedure
with(Data, survival.plot(Procedure, FollowUp, Status, title = "Cancer-Related Death by Surgical Procedure", position = "bottomright", ylim = c(0.5, 1)))

# By histologic grade
with(Data, survival.plot(Grade, FollowUp, Status, title = "Cancer-Related Death by Histologic Grade", position = "bottomright", ylim = c(0.5, 1)))

# By glans invasion
with(Data, survival.plot(Glans, FollowUp, Status, title = "Cancer-Related Death by Glans Invasion", ylim = c(0.5, 1)))

# By coronal sulcus invasion
with(Data, survival.plot(Sulcus, FollowUp, Status, title = "Cancer-Related Death by Coronal Sulcus Invasion", ylim = c(0.5, 1)))

# By foreskin invasion
with(Data, survival.plot(Foreskin, FollowUp, Status, title = "Cancer-Related Death by Foreskin Invasion", ylim = c(0.5, 1)))

# By anatomical level
with(Data, survival.plot(Level, FollowUp, Status, title = "Cancer-Related Death by Anatomical Level", position = "bottomright", ylim = c(0.5, 1)))

# By median tumor thickness
with(Data, survival.plot(Thickness_Median, FollowUp, Status, title = "Cancer-Related Death by Median Tumor Thickness", ylim = c(0.5, 1)))

# By median patient's age
with(Data, survival.plot(Age_Median, FollowUp, Status, title = "Cancer-Related Death by Median Patient's Age", ylim = c(0.5, 1)))

# By urethral invasion
with(Data, survival.plot(Urethra, FollowUp, Status, title = "Cancer-Related Death by Urethral Invasion", ylim = c(0.5, 1)))

# By vascular invasion
with(Data, survival.plot(Vascular, FollowUp, Status, title = "Cancer-Related Death by Vascular Invasion", ylim = c(0.5, 1)))

# By perineural invasion
with(Data, survival.plot(Perineural, FollowUp, Status, title = "Cancer-Related Death by Perineural Invasion", ylim = c(0.5, 1)))

# By pathological T stage
with(Data, survival.plot(pT, FollowUp, Status, title = "Cancer-Related Death by Pathological T Stage", ylim = c(0.5, 1)))

# By tumor relapse
with(Data, survival.plot(Relapse, FollowUp, Status, title = "Cancer-Related Death by Tumor Relapse", position = "bottomright"))

# By local tumor relapse
with(Data, survival.plot(Local, FollowUp, Status, title = "Cancer-Related Death by Local Tumor Relapse"))

# By regional tumor relapse
with(Data, survival.plot(Regional, FollowUp, Status, title = "Cancer-Related Death by Regional Tumor Relapse"))

# By systemic tumor relapse
with(Data, survival.plot(Systemic, FollowUp, Status, title = "Cancer-Related Death by Systemic Tumor Relapse"))

# By clinical N stage positivity
with(Data, survival.plot(cN_Positive, FollowUp, Status, title = "Cancer-Related Death by Clinical N Stage Positivity", ylim = c(0.5, 1)))

# By Solsona et al risk groups
with(Data, survival.plot(Solsona, FollowUp, Status, title = "Cancer-Related Death by Solsona et al Risk Groups", position = "bottomright", ylim = c(0.5, 1)))

# By EAU risk groups
with(Data, survival.plot(EAU, FollowUp, Status, title = "Cancer-Related Death by EAU Risk Groups", ylim = c(0.5, 1)))

# By Hungerhuber et al risk groups
with(Data, survival.plot(Hungerhuber, FollowUp, Status, title = "Cancer-Related Death by Hungerhuber et al Risk Groups", ylim = c(0.5, 1)))

library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
# Plots
par(cex = 1.25)
FN_ROC_Solsona <- with(Data, plot.roc(Final_Nodal, Solsona, main = "Final Nodal Status", lty = 4, col = 4))
FN_ROC_EAU <- with(Data, lines.roc(Final_Nodal, EAU, lty = 2, col = 1))
FN_ROC_Hungerhuber <- with(Data, lines.roc(Final_Nodal, Hungerhuber, lty = 6, col = 2))
legend("bottomright", lty = c(4,2,6), col = c(4,1,2), bty = "n",
legend = c("Solsona et al system", "EAU system", "Hungerhuber et al system"))

# Statistics
ROC_FN_Solsona <- round(with(Data, roc(Final_Nodal ~ Solsona, ci = TRUE)$ci), 2)
MW_FN_Solsona <- format(with(Data, wilcox.test(as.numeric(Solsona) ~ Final_Nodal)$p.value), digits = 2)
ROC_FN_EAU <- round(with(Data, roc(Final_Nodal ~ EAU, ci = TRUE)$ci), 2)
MW_FN_EAU <- format(with(Data, wilcox.test(as.numeric(EAU) ~ Final_Nodal)$p.value), digits = 2)
ROC_FN_Hungerhuber <- round(with(Data, roc(Final_Nodal ~ Hungerhuber, ci = TRUE)$ci), 2)
MW_FN_Hungerhuber <- format(with(Data, wilcox.test(as.numeric(Hungerhuber) ~ Final_Nodal)$p.value), digits = 2)
| Solsona et al system |
0.67 |
0.62, 0.72 |
4.6e-07 |
| EAU system |
0.53 |
0.51, 0.55 |
0.049 |
| Hungerhuber et al system |
0.68 |
0.61, 0.74 |
2.1e-06 |
# Plots
par(cex = 1.25)
FN_ROC_Solsona <- with(Data, plot.roc(DOD, Solsona, main = "Cancer-Related Death", lty = 4, col = 4))
FN_ROC_EAU <- with(Data, lines.roc(DOD, EAU, lty = 2, col = 1))
FN_ROC_Hungerhuber <- with(Data, lines.roc(DOD, Hungerhuber, lty = 6, col = 2))
legend("bottomright", lty = c(4,2,6), col = c(4,1,2), bty = "n",
legend = c("Solsona et al system", "EAU system", "Hungerhuber et al system"))

# Statistics
ROC_DOD_Solsona <- round(with(Data, roc(DOD ~ Solsona, ci = TRUE)$ci), 2)
MW_DOD_Solsona <- format(with(Data, wilcox.test(as.numeric(Solsona) ~ DOD)$p.value), digits = 2)
ROC_DOD_EAU <- round(with(Data, roc(DOD ~ EAU, ci = TRUE)$ci), 2)
MW_DOD_EAU <- format(with(Data, wilcox.test(as.numeric(EAU) ~ DOD)$p.value), digits = 2)
ROC_DOD_Hungerhuber <- round(with(Data, roc(DOD ~ Hungerhuber, ci = TRUE)$ci), 2)
MW_DOD_Hungerhuber <- format(with(Data, wilcox.test(as.numeric(Hungerhuber) ~ DOD)$p.value), digits = 2)
| Solsona et al system |
0.65 |
0.6, 0.69 |
0.00054 |
| EAU system |
0.52 |
0.51, 0.54 |
0.23 |
| Hungerhuber et al system |
0.67 |
0.59, 0.76 |
0.00032 |